home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / t_cvf.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  2.1 KB  |  78 lines
  1. ;$Id: t_cvf.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       T_CVF
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the cutoff value (v) such that:
  11. ;                   Probability(X > v) = p
  12. ;       where X is a random variable from the Student's t distribution
  13. ;       with (df) degrees of freedom.
  14. ;
  15. ; CATEGORY:
  16. ;       Statistics.
  17. ;
  18. ; CALLING SEQUENCE:
  19. ;       Result = t_cvf(P, DF)
  20. ;
  21. ; INPUTS:
  22. ;       P:    A non-negative scalar, in the interval [0.0, 1.0], of type
  23. ;             float or double that specifies the probability of occurance
  24. ;             or success.
  25. ;
  26. ;      DF:    A positive scalar of type integer, float or double that
  27. ;             specifies the degrees of freedom of the Student's t 
  28. ;             distribution.
  29. ;
  30. ; EXAMPLE:
  31. ;       Compute the cutoff value (v) such that Probability(X > v) = 0.025
  32. ;       from the Student's t distribution with (df = 5) degrees of freedom. 
  33. ;       The result should be 2.57058
  34. ;         result = t_cvf(0.025, 5)
  35. ;
  36. ; REFERENCE:
  37. ;       APPLIED STATISTICS (third edition)
  38. ;       J. Neter, W. Wasserman, G.A. Whitmore
  39. ;       ISBN 0-205-10328-6
  40. ;
  41. ; MODIFICATION HISTORY:
  42. ;       Modified by:  GGS, RSI, July 1994
  43. ;                     Minor changes to code. New documentation header.
  44. ;-
  45.  
  46. function t_cvf , a1, df
  47.   
  48.   on_error, 2  ;Return to caller if error occurs.
  49.  
  50.   a = a1 
  51.   if a lt 0. or a gt 1. then message, $
  52.     'p must be in the interval [0.0, 1.0]'
  53.   if (a gt 0.5) then adjust = 1 else begin
  54.     adjust = 0
  55.     a = 1.0 - a 
  56.   endelse
  57.   if a1 eq 0 then return,  1.0e12
  58.   if a1 eq 1 then return, -1.0e12
  59.  
  60.   case 1 of
  61.     df eq 1: up = 100 > (100 * 0.005/a1)
  62.     df eq 2: up = 10  > (10  * 0.005/a1)
  63.     df gt 2 and df le 5:  up = 5 > (5 * 0.005/a1)
  64.     df gt 5 and df le 14: up = 4 > (4 * 0.005/a1)
  65.     else: up = 3 > (3 * 0.005/a1)                
  66.   endcase
  67.  
  68.   while t_pdf(up, df) lt a do begin
  69.     below = up
  70.     up = 2 * up
  71.   endwhile
  72.   
  73.   x = bisect_pdf([a, df], 't_pdf', up, 0)
  74.   if (adjust) then return, -x   $
  75.     else return, x
  76. end
  77.  
  78.